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The sensitivity of trajectories over finite time intervals t to perturbations of the initial conditions can be 
associated with a finite-time Lyapunov exponent A, obtained from the elements Alij of the stability matrix M. 
For globally chaotic dynamics A tends to a unique value (the usual Lyapunov exponent Aoo) as t is sent to 
infinity, but for finite t it depends on the initial conditions of the trajectory and can be considered as a statistical 
quantity. We compute for a particle moving in a random time-dependent potential how the distribution function 
P(A; t) approaches the limiting distribution P{X; oo) — 5{X ~ Aoo). Our method also applies to the tail of the 
distribution, which determines the growth rates of positive moments of Mij. The results are also applicable to 
the problem of wave-function localization in a disordered one-dimensional potential. 
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I. INTRODUCTION 

In this work, we give a uniform description of the complete 
asymptotic statistics of the finite-time Lyapunov exponent for 
a particle moving in a random time-dependent potential. The 
Lyapunov exponent Aoo characterizes the sensitivity of trajec- 
tories to small perturbations of the initial conditions and plays 
a fundamental role in the characterization of systems which 
display deterministic chaos The Lyapunov exponent is 
defined in the joint limits of vanishing initial perturbation and 
infinitely large times. In a hyperbolic Hamiltonian system Aoo 
may be obtained from any non-periodic trajectory, because for 
arbitrarily long times the trajectories uniformly explore the 
complete phase space. 

A widely studied generalization of Aoo is the finite-time 
Lyapunov exponent [|,|,|,|,||0|g|ig|ll],g|l3[ 
p^, [Is, 16, 17 1, which is defined for finite stretches (time in- 



terval t) of trajectories (generalizations to finite perturbations 
also exist JT^). The sensitivity of the dynamics to initial per- 
turbations is given by the stability matrix map M, which is the 
linearization of the map of initial coordinates to final coordi- 
nates. In terms of elements Mij of M, the (complex) finite- 
time Lyapunov exponent may then be defined as 



A = -InMj, 
t 



(1) 



In contrast to Aoo, A is not a unique number independent 
of the initial conditions, but a fluctuating quantity with a dis- 
tribution function P(A; t) (defined by uniformly sampling all 
initial conditions in phase space). This distribution function 
determines, e. g., the generalized entropy and dimension spec- 
tra of dynamical systems and more practically the weak- 
localization correction to the conductance [ p^ and the shot- 
noise suppression [^plt| inmesoscopic systems. Finite-time 
Lyapunov exponents also determine the wavefront stability of 
acoustic and electromagnetic wave propagation through a ran- 
dom medium, in the ray-acoustics/ray-optics regime of short 



wave lengths (for a recent application see Refs. [ [22| , |23[]). 
Moreover, they have shifted into the focus of attention due 
to recent advances in the understanding of the role of the 
Lyapunov exponents for quantum-chaotic wave propagation 
Ig, H H, I?! ^ ||,|0||: It has been observed that under 
certain conditions the Lyapunov exponent can be extracted 
from the decay of the overlap of two wavefunctions which 
are propagated by two slightly different Hamiltonians (the so- 
called Loschmidt echo). Since the overlap is studied as a func- 
tion of time, the distribution of the finite-time Lyapunov expo- 
nent is directly relevant for these investigations. This extends 
also to related semiclassical time scales, like to the Ehrenfest 
time ~ (log?i)/A, which is a semiclassical estimate of the 
diffraction time of wave packets due to the chaotic classical 
dynamics. 

In the limit of infinite time t the distribution function 
P(A; <) in a completely chaotic phase space tends to the limit- 
ing form P( A; cx)) ~ ^(A — Aoo)- For large but finite t the bulk 
of the distribution function can be approximated by a Gaus- 
sian centered around Aoo, with the width vanishing cx t^^/^ 
as t oo. However, many of the properties determined by 
P(A; t) (like the generalized entropy and dimension spectra) 
cannot be calculated from the Gaussian bulk of the distribu- 
tion function [|l]] . 

In this paper we investigate for a particle moving in a one- 
dimensional random time-dependent potential how P(A; t) 
approaches the limiting distribution function P(A; oo) = 
S{X — Aoo) for large times. Our approach uniformly applies 
both to the bulk as well as to the far tail A ^ Aoo of the 
distribution function. We find that the cumulant-generating 
function of P(A; t), 

viO = ln(exp(etA)) = ^((A"))iif , (2) 

11=1 

(where the average (•) is over initial conditions and ((•)) de- 
notes the cumulants), takes the asymptotic form 

viO^KOt/tc + o{f), (3) 

with ^(^) a universal function (within the statistical model) 
and 
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tc 



A^oV'^^ 



(4) 
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a system-specific time-scale which can be determined from 
the infinite-time-Lyapunov exponent and the constant /i*^^' = 
d/i/(i^|{=o (by definition, dri/d^\^^() — Xoct)- The function 
/i(^) is given by the leading eigenvalue of a second-order dif- 
ferential equation in which ^ appears as a parameter This 
eigenvalue can be calculated perturbatively in ^, which gen- 
erates the cumulants of A. The values of /i at integer ^ de- 
termine the asymptotic growth rates (l/t) ln(Af|^) = /i(C)/^c 
of positive moments of elements of the stability matrix. We 
find that these values are given by the leading eigenvalue of 
finite-dimensional matrices. 

A random time-dependent potential is often considered as 
a statistical model for the ergodic properties of hyperbolic 
chaotic motion, in the spirit of the early work of Chirikov pi] ] . 
The time dependence of the potential may be considered to 
mimic the dependence of the potential in the eigentime along 
the trajectory. In the context of finite-time Lyapunov expo- 
nents there have been indications that a statistical description 
is usually valid for the chaotic background of its distribution 
JT^], while system-specific deviations may exists in some ex- 
ceptional cases even in the bulk of the distribution function 
Jl^]. While the statistical model considered in this work is 
tailored to a specific class of Hamiltonian systems, it can be 
modified straightforwardly to other classes of chaotic systems 
(this is briefly described at the end of this paper). 

The problem of finite-time Lyapunov exponents in the ran- 
dom time-dependent potential is equivalent to the problem 
of wave-localization in a random one-dimensional potential 
p3| , p4| , ^ 36 1, because the equations of motion for 
the matrix elements Mij are formally equivalent to the corre- 
sponding Schrodinger equation JlC| , p7| ]. Indeed, the Fokker- 
Planck equation employed in this work is based on the phase 
formalism described, e.g., in Ref. [Bsl 59l 40|. Hence, the 



p), and the Hamiltonian given by 



asymptotic statistics of the finite-time Lyapunov exponent 
presented in this work directly is of interest and can be trans- 
ferred to this field of research. A number of additional areas 
of application of our method come into scope if one considers 
the vast arena of problems which can be analyzed by products 
of random matrices, since the finite-time Lyapunov exponents 
are a valuable way to characterize the eigenvalues of these 
products [[To|]. 

The plan of this paper is as follows: In Sec. ^ we for- 
mulate the problem of finite-time Lyapunov exponents in the 
one-dimensional random time-dependent potential. In Sec. |l| 
we show how the cumulant-generating function can be related 
to the parameterized eigenvalue of a second-order differential 
equation, and that the cumulants can be calculated systemati- 
cally. Positive moments of Mij are calculated in Sec. |^ We 
close the paper with discussion and conclusions in Sec^Vj 



II. FORMULATION OF THE PROBLEM 

A. Statistical model 

Let us consider a time-dependent Hamiltonian system with 
one degree of freedom (canonically conjugated coordinates x. 



(5) 



Here V{x, t) is a time-dependent potential and m is a mass. 
We also allow for an additional static potential with curvature 
V2 acting in the background of the random potential (this po- 
tential is repulsive for V2 < and attractive for V2 > 0). 

We introduce the map Tt which propagates initial condi- 
tions {xi,pi) over a time interval t to the final coordinates 
{xf,pf) — Tt{xi,pi). The stability matrix A/ is the lineariza- 
tion of the map Tt and describes the sensitivity of the final 
coordinates to a small perturbation of the initial conditions. 



M 



9{xf,Pf) 
dixt,Pi) 



Mil M12 
M21 M22 



(6) 



Area preservation of the dynamics in phase space entails the 
property det M = 1 of the stability matrix. 

We are interested in the evolution of the stability matrix 
with given initial conditions and increasing time interval t. 
According to Hamilton's equations of motion the stability ma- 
trix fulfills the differential equation 



dM 
~dt 



= KM, K = 



m ^ 
V 



where the function v{t) in the matrix K is given by 

d'^V 



-V2 



dx"^ 



(7) 



(8) 



{x,p) = {xf,pi) 



This differential equation is supplemented by the initial con- 
ditions 



M(0) =diag(l,l), 



(9) 



corresponding to the identification of the initial and final co- 
ordinate systems for < = 0. 

In order to study the statistical behavior of the stability ma- 
trix we now assume that v{t) is a randomly fluctuation func- 
tion equivalent to Gaussian random (5-correlated noise, 

{v{t)) = -V2, {viti)vit2)) = 2D6{ti - t2). (10) 

The condition of a vanishing mean of the time-dependent part 
of V corresponds to the observation that the incidence of pos- 
itive and negative curvature of the potential landscape along a 
typical chaotic trajectory should be identical. The (5-function 
correlations are valid if the correlation time of the fluctuations 
is smaller than the mean free transport time in the random po- 
tential. The constant D (similar to a diffusion constant, but 
not identical with conventional diffusion constants of motion 
in phase space) can be related to the strength of the temporal 
fluctuations of the potential V{x, t). However, both D as well 
as the mass m can be eliminated from the subsequent analysis 
by rescaling quantities in the following way: 

t = tct', v^{D/m)tcv\ V2^{D/m)tcV2 
M12 = itc/m)M[2, M21 = (m/QM^i, 
Mii=M'ii, M22 = M^2- (11) 
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Here we defined the characteristic time scale 



(12) 



[In the course of our analysis we will see that this time scale 
also can be found from Eg . (H ).] The rescaled (primed) quan- 
tities fulfill Eqs. (0), (|), with D = m = 1. Also note 
that the rescaling leaves the property det M — 1 invariant. 



B. Relation to one-dimensional localization 

The set of linear first-order differential equations (Q) can be 
decoupled by converting them into second-order differential 
equations. It is useful to note (as mentioned in the introduc- 
tion) that the equations for the elements Mn and M12 are 
equivalent to the Schrodinger equation, at energy E = V2/ m, 
of a particle of mass /2 in a one-dimensional random po- 
tential (v + V2)/m (of vanishing mean), with t playing the 
role of the spatial coordinate. 



d^Mu 



V 

-Mil, 
m 



V 

-M12, 
m 



j,n — ^'^llj ,,n — ^'-iizi (13) 

dt^ m dt'' m 

while the other matrix elements are directly related to them by 



M21 = m 



dMii 
dt 



M22 = m- 



dM^ 



12 



dt 



(14) 



The problem of finite-time Lyapunov exponents hence is 
closely related to the problem of one-dimensional localiza- 
tion in a random potential, in which the Lyapunov exponent 
corresponds to the inverse decay length of the wave function. 



III. CUMULANTS OF THE FINITE-TIME LYAPUNOV 
EXPONENT 



We now solve the problem of finding the probability distri- 
bution function of matrix elements My within the statistical 
model of chaotic dynamics, defined by the evolution equation 
(0) for M, with initial condition (Q), and the statistical proper- 
ties ( [l0| ) of the random function v. For the sake of definiteness 
we will consider in this section the statistics of the upper diag- 
onal element Mn. The results directly carry over to the other 
elements of M, as is discussed in Sec. 



Note that the evolution equation of z decouples from u and 
can be interpreted as a Langevin equation. Hence the distribu- 
tion P{z\ t') can be calculated from a Fokker-Planck equation, 
which was considered before in the context of wavefunction 
locahzation ||3|,|^, 

dt'P{z-t') = C,P{z-t'), (17a) 
■ = d,{z^ + Vi + d,) ■. (17b) 

For large t' the distribution function P(z;t') approaches the 
stationary solution [|8l |9|, |o| 



PU^) = N / dyK{y,z) 



N 



(18a) 



(18b) 



7T-^[Ai'{-V^)+Bi\-V^)]-\ (18c) 



Here Ai and Bi are Airy functions. The normalization con- 
stant is directly related to the integrated density of states in 
the localization problem @ ^ |oll. For ^3' = 0' ^ = 
35/62-i/37r-i/2/[r(l/6)]. Because du/dt = z/U it is cleai- 
| |4C| ] that the infinite-time Lyapunov exponent can be obtained 
from Aoo — j z) /tc\ this relation will be demonstrated ex- 
plicitely in Sec. [11 C . 

The Fokker-Planck equation for the joint distribution func- 
tion z\ t') is given by 



dvP = -zduP + C,P. 



(19) 



This Fokker-Planck equation with V2 = has been derived 
in Ref. [ [l9| ] for the autonomous chaotic scattering of a parti- 
cle from a dilute collection of scatterers (with more than one 
degree of freedom). 

The joint distribution function P{u, z; t') does not approach 
a stationary limit because u runs away to infinitely large val- 
ues. In order to analyze the behavior of the distribution func- 
tion P{u, z\ t') for large times we convert the Fokker-Planck 
equation ( p9| ) into an eigenvalue problem which discriminates 
between the different time scales involved in this evolution. 
For this purpose, we introduce into Eq. ( [l^ ) the ansatz 



P{u,z;t') = 



/+ioo 
-too 27ri f-' 



u)U^,z). (20) 



A. Cumulant-generating function as an eigenvalue 



We introduce the quantities 

u = InM^i, 



M|i 



(15) 



where the relation u = At to the finite-time Lyapunov expo- 
nent A is established by Eq. (|l|) [note that Afn = M[i in the 
rescaling Eq. ([TTI)]. According to Eqs. (^ and ([TTl), u and z 
fulfill the differential equations 



du 
dp 



dz 
dt' 



v' - z\ 



(16) 



(The integration contour along the imaginary axis corresponds 
to an inverse Laplace transformation.) It follows that the func- 
tions /„ fulfill the differential equation 



Mn/n(C, Z) = {(Z + C^)f„iC, z), 



(21a) 



in which /i„ and ^ appear as parameters. However, in order to 
obtain a meaningful probability distribution function (20) we 
have to impose boundary conditions on fn{^, z) at z ^ ±00. 
It is convenient to express these boundary conditions by the 
requirement 



dz fn{^, z)z < 00. 



(21b) 
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Here V denotes the principal value w ith re spect to the inte- 
gration boundaries at ±00. Condition ( 21b ) follows from the 
behavior z k, {t' — t'^)~^ of the solution of the differential 
equation ( [1^ close to times t' « t'^ where |z| ^ cx3 (and 
hence v' can be ignored). In practical terms, the condition 
( 21t ) guarantees that the drift of u remains finite for all times. 

Eq s. ( pT| ) form an eigenvalue problem, since condition 
(2 It) only can be fulfilled for a discrete set of numbers /i„ — 
note that these eigenvalues depend on the parameter ^. In the 
limit of large t' only the largest eigenvalue ^o(0 = l^iO 
is relevant, because the other eigenvalues give rise to expo- 
nentially smaller contributions. This eigenvalue vanishes as 
f ^ 0, /. e., — 0, because the stationary distribution of 
z, Eq. (p^, must be recovered for large times from Eq. ( pO[ ) 
by integrating out u. 

The moments of u are given by 



+^oo ^ 



du exp(/if' — £,u)f{^)u^' 



(22) 



where the coefficients /(^) = J^^dzfo{£,, z) are deter- 
mined, in principle, by the initial condition for P{u, z; t') at 
t' = 0. From Eq. ( p2[ ) we obtain the moment-generating func- 
tion 



X(0 = {eM^u)) = exp(/i(e)t/tc)/(C), 



(23) 



where we re-introduced the original time variable t — tct' by 
Eq. ([TTI). The cumulant-generating function (^ hence takes 
the form of Eq. (^, including the corrections of order t^. 



viO = ^^x{0^f^{Ot /tc + In fiO- 



(24) 



The cumulants ((A")) of the finite-time Lyapunov exponent 
are obtained by expanding the generating function rj in powers 
of ^, see Eq. (^. In terms of the coefficients of the Taylor 
expansion 



in) 



(25) 



n=l 



[which starts with the linear term in ^ because /i(0) = 0], 
according to Eqs. (^ and ( |24| ) the nth cumulant of A is then 
given by 



((A"» 



o{t- 



(26) 



This equation means that within the statistical model the cu- 
mulants are universal quantities in the leading order in t, in the 
sense that the initial conditions P{z, u; 0) only enter the next- 
order corrections. The only system-specific parameters which 
enter the cumulants are the time scale tc and the (rescaled) 
strength V2 of the static potential. Note that ratios of cumu- 
lants are even independent of the time scale tc (and hence of 
the parameters D and m of the statistical model). 

The form (^ of tc follows from Eq. (^) when tc is ex- 
pressed in terms of the infinite-time Lyapunov exponent Aoo 
with help of the definition 



- „(i) 



Aoo = lim (A) = /i 

t — >oo 



/tc 



(27) 



In terms of the bare quantities of the statistical model, 

= ^(l)i?l/3^-2/3_ (28) 

In the next two sections we obtain general expressions for 
the expansion coefficients /i'") and calculate explicitely the 
proportionality factor /x'^^ — d/i/d^j^^o in (p8|), as well as 
the first few coefficients /i^^-', . . . , which determine, re- 
spectively, the variance and the leading non-Gaussian correc- 
tions (higher cumulants) of the fluctuations of the finite-time 
Lyapunov exponent around its limiting value Aoo ■ 

B. Recursion relations for the cumulants 

We now show how the cumulants can be calculated from 



Eq. (26) by recursively solving a hierarchy of equations for 
coefficients fi^"'' in the Taylor expansion of /i(^), Eq. (^5|). 

In analogy to Eq. ( p5| ) let us also expand the function 
fo{£.,z) in powers of 



(29) 



With Eqs. ( |25| ) and ( |29| ) the eigenvalue problem ( |21[ ) can now 
be written order by order in powers of For n = we 
recover the stationary variant ([Tt]) of the Fokker-Planck equa- 
tion 



(30) 



which is solved by the stationary solution fQ^\z) = Pstat(-z), 
Eq. ( 18 ). For n > 1 the differential equations are of the form 



(z) = -zft^'\z)+yu(yt-'\z) 



(31) 



Let us assume that we have solved the hierarchy of equa- 
tions up to order n — 1. In the next order n both the unknown 
quantities /q""* as well as /i*^"^ appear The unknowns can be 
separated by integrating the differential equation ( ^1] ) over z 
from —00 to 00: The i ntegrated left-hand side vanishes be- 
cause of condition ( [21b| ) of the eigenvalue problem. The inte- 
grated right-hand side can be rearranged to give /i^"-'. 



,(") 



/oo «-l 
dz[zf^-'\z)-Y^^,(y^-'\z)], (32a) 
1=1 



which only involves quantities up to order n — 1. Subse- 
quently, /x*^"' can be inserted into Eq. (pi}). The function 



dy / dxK{y, 



x[-xft-'\x)+Y,^^^^^ft^'\x)] (32b) 
1=1 

[with the kernel K{y, z) defined in Eq. ([T^] is then obtained 
by solving the resulting inhomogeneous differential equation 
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with help of the partial solution flf\z) of its homogeneous 
counterpart, Eq. (BQ). This inhomogeneous part of the func- 
tions /g"' (z) is fixed by the requirement that /q"-* (z) is nor- 
malized to 1. Adding the homogeneous solution to f^^\z) in 
any order gives rise to additional terms in all higher orders, 
but these combine in such a way that they drop out of the 
calculation of the c oeffi cients /i'"-*, which hence are uniquely 
determined by Eq. (32a). 



The recursion relations (32 ) can be iterated to calculate suc- 



cessively all cumulants of A. 




C. Explicit expressions and numerical values 



According to Eq. (Pq), the two numbers ^i^^^ and /i^^^ de- 
termine mean and variance of the distribution function of A, 
which then is approximated by a Gaussian. The coefficient 
jjS^^ has been obtained in Ref. from the Fokker-Planck 
equation ( p^ for arbitrary V2. For the special case V2 — 0, 
the two coefficients /i*^^^ and pS'^'^ have been obtained in Ref. 
JT^ ] from the Fokker-Planck equation (19). However, the de- 
viations from the Gaussian distribution function are not at all 
negligible for many chaotic systems, which is most clearly 
displayed in their generalized dimension and e ntropy spectra 
[jl|]. As we have seen in the previous subsection |III B| , our ap- 
proach of reduction to the eigenvalue problem (|21[) allows to 
analyze the non-Gaussian deviations by the higher cumulants 
of A. [In next section IV, we show that one can even obtain 
from our analysis the positive moments of Afn, which are de- 
termined by the far tail A ^ Aoo of P(A; t), while the bulk of 
the distribution is essentially irrelevant for these moments.] 

,(2) 



Explicit expressions for the first few coefficients 
/i^'^), and /x^^^ result from Eq. (32a), 



,(1) 



dz{z-ix(^'>)f^^\z), 



1^' 

(33a) 
(33b) 



,(3) 



,(4) 



(33d) 



where flf\z) = fstat(^) is given by the stationary distribu- 
tion function of z, Eq. (|l8|), while the other functions follow 



=1. 
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FIG. 1: Coefficient /^'^^ of the first cumulant, and the ratios 
n!//'"'//^'^' for the coefficients of the second, third, and fourth cu- 
mulant [cf. Eq. ([2^], as a function of the strength Vi of the static 
background potential. 
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84^/^ 




5 


6 


7 




2(14 + 3^)^/^ 


(252 + 24^)^/^ 


2(63 + 15^10)^''^ 



TABLE I: First eight coefficients n|u'"' of the cumulants of finite- 
time Lyapunov exponents [cf. Eq. (Eq)], in absence of the static back- 
ground potential (V2 = 0). 



TABLE II: Exponential growth rates /i(^) of the first few moments 
(A/jj) [cf. Eq. (^^], in absence of the static background potential 
(Vi = 0). 



from Eq. (32b), 



dydxK{y,z){fi^^'> -x)f^°\x), (34a) 

z>y>x 
z>y>x 



(34b) 



dydxK{y,z)[{^l^'^ -x)fl,'\x) 

z>y>x 

(34c) 



The coefficient /i^^^ is then given by 



(35) 



where N is given in Eq. (|18|), while the cumulants for n > 2 
can be obtained quickly by numerical integration of 2ri-fold 
integrals. The effort of integration can be greatly reduced 
down to the expense equivalent to a single integral, because 
the integrand factorizes. An efficient recursive scheme is de- 
scribed in the Appendix. In Fig. |l] we plot the coefficient /i'^^ 
and the ratios for n — 2,3,4 as a function of 

V2 . The non-Gaussian corrections are largest around V2 = 0, 
while they become irrelevant for large negative or positive val- 
ues of ¥2- 

For strong confinement, V2 ^ 1, the coefficients 
n!/i'^"^//i^^) Sin + ^2n^ with Smn the Kronecker symbol, 
and the Gaussian approximation 



(36) 



becomes valid. [In the context of wave localization, this cor- 
responds to the well -known limit of a large Fermi energy 
E - (cf. Sec. flB]).] 



Analytical results can be found in the case V2 = for the 
first two coefficients. 



^(1) = 



(3/2) 



1/3, 



r(i/6) 



18 



-N - 



7 3 3 3 

' '6' 2' 2' 4 



(37a) 
(37b) 



where iV(F2' = 0) = 35/^2-i/37r-i/V[r(l/6)], while gi^a is 
a generalized hypergeometric function. Incidentally, the nu- 
merical value given for /i*^^^ in Ref. [|l9| is wrong, but the 



analytic e xpres sion given in that paper is equivalent to Eqs. 
( 33b ) and (37b). In Tab. | we tabulate the numerical values of 
the first eight coefficients for = 0. 



IV. POSITIVE MOMENTS 

A. Formally exact expressions 

In view of Eqs. (|l5]) and ( p3[ ) we find that the exponential 
growth rates of the positive moments of Mn are given by the 
eigenvalue /x(^) of Eq. ( ^If ) at integer values of ^: 



dt 



tc 



(38) 



As we will now show, for integer values of ^ the eigenvalue 
problem (|T]) can be reduced to a matrix eigenvalue problem 
of finite dimension. For the first few moments the leading 
eigenvalue can be calculated explicitely, while for larger val- 
ues it is formally given by the largest root of the corresponding 
characteristic polynomial. 

In order to obtain a solution of the differential equation 
(21a), we write 



dyK{y,z) 



9i± 

9{yf 



(39) 



[with the kernel K{y, z) defined in Eq. (|18[)], and obtain for g 
the differential equation 



{^l-^,z)g = -{z^ + Vi)^,g + ^lg 



(40) 



(a triconfluent Heun's equation with singularity at 1/z — 0). 
We introduce into this equation the polynomial ansatz 



(41) 



Power matching results in the following recursion relation 



{n + 3)c„+3] 

(42a) 



for the coefficients c„, with initial conditions 

= 1, Cj_i = C^-2 = M^/2. 



(42b) 
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FIG. 2: Growth rates /i(^) of the moments {M^i) [cf. Eq. @] in 
absence of the static background potential (V2 = 0), obtained as the 
largest eigenvalue of the matrix (H) (full circles). Also shown is the 
real part of the subleading eigenvalue of this matrix (open circles). 



For integer ^ this recursion relation terminates. We obtain 
functions co(/x), ci(/i), and C2{^) and an additional condition 
from the term in Eq. ( p^ which is constant in z. 



p^{n) = ^co + V2C1 - 2c2 = 0, 



(43) 



where ^^(/i) is a polynomial of degree f + 1. 

The polynomial (/i) can also be interpreted as the char- 
acteristic polynomial of the + 1) x + 1) -dimensional 
matrix 



{ -Vi 1-2 









f -2V^ 2-3 
C-1 -Wi 3-4 
-WJ, 

C-3 



(1-0^2 (e-i)e 







(44) 



which is the matrix representation of the eigenvalue problem 
( pll ) in the space of the monomial expansion of g{z). 

The exponential growth rate of the ^th moment is 
given by the largest root of (/i), or equivalentl y by the 
largest eigenvalue of the matrix (H). In subsection [VB we 
will see for the examples ^ = 1,2 that the other roots show up 
in the transient behavior of the moments. 

First we present results in absence of the static background 
potential, V2 = 0. The values for the first few moments are 
given in Tab. 0. Figure || shows the growth rates and the real 
part of the subleading eigenvalue for values of ^ up to 80. 
A log-normal statistics of Afn (corresponding to a Gaussian 
statistics of the finite-time Lyapunov exponents) would result 
in the quadratic dependence Eq. (|6|) of /i(^) on ^, while the 
plot shows a weaker (approximately linear) dependence for 



large ^. This results from the influence of the terms /^("V" 
for n > 3 in the complete Taylor expansion of /x, Eq. (g^). 
Further note that the subleading eigenvalue stays at a finite 
distance to the leading eigenvalue (indeed, their distance in- 
creases with increasing ^), as we have assumed before in 
restricting the attention to the leading eigenvalue //q of the 
eigenvalue problem (pi)). 

For finite V2, the growth rate of the first moment 



Re/^(1) 



vanishes in the case of confinement, V2 > 0. This will be 



confirmed by the direct computation in Sec. 
rate of the second moment is given by 



IV B 



(45) 

11 be 
The growth 



fi{2) = 2^/3 + ;/l + 16^2' 

2^/^ ( 1 - a/1 + 16y2'727 



1/3 



1/3 



(46) 



[with the roots taken such that /i(2) is real]. We plotted the 
real parts of the leading and subleading growth rates [eigen- 
values of matrix (p4|)] for the first four moments in Fig. 



B. Direct computation of the first and second moment 

In order to illustrate the results for the growth rates of the 
moments (Af^j^) we compare the results for ^ = 1 and ^ = 2 
to the exact results for all times (including the transient be- 
havior). A formal solution of the differential equation ( 13 ) in 
terms of a series in the disorder potential is obtained by inte- 
grating Eq. (|l3]) twice, under observation of the initial condi- 
tions Mil ~ 1^ dMii/dt — for t — 0, and iterating the 
resulting integral relation 



Mii{t) = 1 



* dti dsi^^Miiisi) 











1+ f dti{t-tif-^Mii{ti] 
Jo m 



The formal solution is of the form 



Mn(to) = l + En 



■11=1 k=i 



dtk (t 



tk) 



v{tk) 



(47) 



(48) 







where we introduced = t for notational convenience. 

For the first moment we can average Eq. (^8|) directly. Be- 
cause of the factors {tk-i — tk) and the time ordering, the 
random function v never appears instantaneously in second or 
higher order in any of the integrals. Hence we can replace v by 
its average, given in Eq. (|l^. Consequently, the first moment 
is given by 



(Mil) =cos[(t/t,)VW 
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For V2 = the first moment is constant and given by its ini- 
tial value, (A/11) — 1- This means that negative deviations 
Mil ^ 0, corresponding to inverse hyperbolic motion, can- 
cel precisely the positive deviations AIn ^ of hyperbolic 
motion. For negative V2 the first moment grows, while for 
positive V2 it oscillates and stays of order unity. In the decom- 
position of the cosine into the two exponentials, we identify 
in the exponents the two roots ± \/—V^ of the characteristic 
polynomial p^=i(/z) = /z^ + V2 of the matrix ( p4| ) with ^ = 1. 
For negative V2, the subleading exponent hence governs the 
transient behavior of the first moment. 

For the second moment let us restrict for simplicity to the 
case V2 = 0. We group the functions v in the two factors 
of Mil in pairs and then invoke the delta-correlations of Eq. 
Jicl). Performing the time-ordered integrals we obtain 



(A'^n> = i + E(2Acrn 



n=l 



k=l 



dtk {tk-1 — tkY 



(50) 



The asymptotic growth rate of the second moment is given 
by the leading root /i(2) = 2^/'^ of the characteristic pol 
nomial P2(m) = ^fJ.^ — 2, which is in accordance to Tab, 
The second and third exponent are the other two roots of this 
polynomial. 



C. Equivalence of matrix elements 

So far we mainly studied the statistics of the upper diagonal 
element Mu of the stability matrix M. At this point now 
we can discuss how the results for the cumulant-generating 
function and the positive moments can be transferred to the 
other elements of M. 

The differential equation ( |l4| ) for A/22 can be integrated 
similarly as the one for Mu, from which we obtain analo- 
gously to Eq. (Ssh the formal solution 



M2 



n=l k=l 



dtk {tk — tk+i) 



v{tk) 



(51) 



Here we defined in each term of order n that i„+i = 0. It 
follows by direct computation that the first two moments of 
M22 are identical to those of Mu, 



{M22) = (Mil), (MI2) = (Mi\). 



(52) 



These explicit results already suggest that the statistics of the 
two diagonal matrix elements is the same. Indeed, the trans- 
formation ifc = t—tn+i-k, v{t—t) = ^(t) brings Eq. (51 ) into 
the form of Eq. (^sj) and leaves the properties of the Gaussian 
noise (llOb invariant. Hence even the transient behavior of the 



diagonal elements is completely identical, for arbitrary values 
of VI 

The results for the cumulant-generating function r/(^) 
(hence also the growth rates of the moments, but not the tran- 




FIG. 3: Growth rates ^i(^) the moments (M^) [cf. Eq. (|3^], for 
m = 1,2, 3, 4, as a function of the strength V2 of the static back- 
ground potential. Also shown (dashed lines) is the real part of the 
subleading growth rates [subleading eigenvalues of matrix (M)]. 
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sient behavior) can also be transferred to the offdiagonal ma- 
trix elements of M: The element A/12 fulfills the same dif- 
ferential equation as Mn, see Eq. (|l3|), while M21 fulfills the 
same differential equation as A/22- The initial conditions of 
the offdiagonal matrix elements differ from those of the di- 
agonal elements. However, according to Eq. ( |24| ) this only 
affects the function /(^) in the subleading corrections of the 
cumulant-generating function [which, for the example of the 
second moment, results in factors in front of the exponential 
functions which are different than in Eq. (^0|)]. 

Let us add that from Eqs. and we find for V2 = 
the cross-correlator 



{M,,M22) = \ + \{Ml,). 



(53) 



As a consequence, for V2 ~ the trace tr A/ = A/n + A/22 
of the stability matrix has the following first two moments 

(trA/) = 2, (54a) 
((tr A/)2> = 1 + e''(2)t/t. ^ 2 Ree-("l)''''^(^)*/*^ 

(54b) 



V. DISCUSSION 

In this work we presented a uniform approach to the asymp- 
totic statistics of finite-time Lyapunov exponents, for the 
model (described in Sec. ||) of a particle moving in a random 
time-dependent potential. The cumulant-generating function 
r7(^) was found to be directly proportional to the eigenvalue 
of a parameterized differential equation, defined by Eqs. 
(|2l|). This facilitated an effective analysis of the statistics, in- 
cluding the non-Gaussian deviations of the distribution func- 
tion. These deviations are especially important for the positive 
moments of the elements of the stability matrix, since their 
growth rate cannot be predicted by the Gaussian approxima- 
tion Eq. (36). 

We limited our attention to the case of time-dependent 
Hamiltonian systems with a single degree of freedom and a 
Hamiltonian (^) which is of the special type of kinetic en- 
ergy plus potential energy, with time-dependence only in the 
potential energy. This case is of particular interest because 
of its direct applicability to specific dynamical systems as in 
the random wave-propagation problem of Refs. ^3[], and 
because of its applicability to one-dimensional wave localiza- 
tion. For the Hamiltonian the matrix K in the differential 
equation (^ is purely off-diagonal, with fluctuations only in 
the lower-left element. For Hamiltonians which do not sepa- 
rate into kinetic and potential energy, the differential equation 
P) for A/ involves the matrix K in the more general form 



K 



( Kn 


Ki2^ 




V K21 


K22 ) 





d-'H 
dxdp 



dxdp 



(55) 



A generalized statistical model now arises by introducing 
noise into all of the matrix elements of K. (One may also al- 
low for correlations between the different matrix elements or 



for finite correlation times by introducing auxiliary variables 
for the noise in the standard way.) 

Let us point out two particular cases for which a statis- 
tical description promises to result in direct applications to 
physical situations of interest. One case is more relevant to 
wave-function localization while the other is more relevant for 
chaotic dynamics. 

a) The diagonal elements Kn = —K22 = still vanish 
identically, but both off-diagonal elements K12 and K21 fluc- 
tuate with a vanishing mean. This situation appears to be re- 
lated to the band-center case of one-dimensional localization 
in the Anderson model [ pT[ Q (where space is discretized on 
the lattice), since at the band-center the effective mass of the 
particle diverges (and hence the mean of K12 vanishes). 

b) Chaotic dynamics with an isotropic phase space may 
be modeled by independent fluctuations of all four matrix 
elements Kij with identical amplitude and vanishing mean. 
Hamiltonian dynamics gives rise to the further constraint 
Kii ~ —K22- Isotropic dynamics arises in typical chaotic 
maps (some maps, like the Baker map or the cat map, how- 
ever, are not isotropic — the directions of stable and unsta- 
ble manifolds are known by construction). Good candidates 
are the Poincare surface of section of autonomous systems 
with two degrees of freedom, in which the motion in four- 
dimensional phase space is restricted to three-dimensional 
manifolds of constant energy and the coordinate along the 
flow field is taken as a time. 

It would be interesting to compare the outcome of an anal- 
ysis of model b) with the findings in the literature [15 |l6|] 
which indicate a certain degree of robustness (if not univer- 
sality) of the distribution of finite-time Lyapunov exponents 
for typical chaotic systems. 
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APPENDIX: INTEGRALS FOR THE HIGHER CUMULANTS 

The cumulants of order n result from the recursion rela- 
tions Eq. ( ^2[ ) in the form of 2n-fold integrals. Usually, the 
numerical evaluation of such integrals for large n is very time- 
consuming, since the number of points on a grid covering 
the integration domain with lattice constant (1/A^), N ^ 1, 
grows rapidly with n as A^^". However, presently the in- 
tegrand factorizes and the expense of the integration can be 
reduced from exponential to algebraic n-dependence ^ nN. 
The principle can be demonstrated for the example of the two- 
fold integral 

1 J -za 

(56) 
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where g is an arbitrary function and /^'^^ may itself be a multi- 
dimensional integral. 

We introduce an index m which denotes that the argument 
of a function is taken at the mth lattice point on the appropriate 

(n) 

axis of the grid. The initial values of /,„ ' at to = (the lower 
integration boundary) are zero. We now can write recursively, 
by incrementally increasing the integration variables, 



r(2) 
'-m+1 

'm+1 



9ri 



±/(2) 



9m 



1 

N 



a 



r(2) 



N 



(57a) 
(57b) 



Moreover, when I^^^ itself is a multi-dimensional integral of 
type /(^\ its current value can be obtained recursively in the 
same way as the value of Since each additional integral 
will gi ve ri se to only one additional equation [similar either 
to Eq. (57a) or to Eq. ( |57b| )], the number of operations grows 
linearly with n, as advertised above. [The recursion relations 
( |5^ ) have the additional advantage for the present problem that 
they avoid over- and underflow in the evaluation of the kernel 
Kiy, z) = exp(j/3/3 + V^y - z^/S - V^z).] 
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